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Abstract 

We report the diffusion coefficient and viscosity of popular rigid water models: Two non- 
polarizablc ones (SPC/E with 3 sites, and TIP4P/2005 with 4 sites) and a polarizable one (Dang- 
Chang, 4 sites) . We exploit the dependence of the diffusion coefficient on the system size [Yeh and 
Hummer, J. Phys. Chem. B 108, 15873 (2004)] to obtain the size-independent value. This also 
provides an estimate of the viscosity of all water models, which we compare to the Green- Kubo 
result. In all cases, a good agreement is found. The TIP4P/2005 model is in better agreement 
with the experimental data for both diffusion and viscosity. The SPC/E and Dang-Chang water 
overestimate the diffusion coefficient and underestimate the viscosity. 



I. INTRODUCTION 



The primary importance of water in biological, environmental or technical processes ex- 
plains the enormous efforts devoted to its modelling on the microscopic scale. Many classical 
force fields for molecular simulations have been introduced to capture different aspects of its 
complex behaviour, including a number of " anomalies" . These models vary in the number 
of considered interaction sites, their flexibility or rigidity and their treatment or neglect of 
polarizability. Reviews of water models^"^- have addressed their relative merits to model the 
thermodynamic (e.g. phase diagram or heat capacities), structural and dynamic properties 
of water. The latter usually include the diffusion coefficient and the shear viscosity of the 
model, which are key quantities in determining the transport of water and solutions at rest 
or under hydrodynamic flows. In their recent review of the properties of rigid water models, 
Vega and Abascal underlined that dynamic properties are very sensitive to the details of 
the water model and that "the diffusion coefficient has not received the attention it should 
deserve as a target property" in developing potential modeld^. 

The computation of diffusion coefficients from molecular dynamics simulation must be 
done carefully, as the result may depend on the simulation parameters (cut-off radius or 
method to compute electrostatic interactions), but also exhibits a systematic size-dependence 
arising from the screening of hydrodynamic flows under periodic boundary conditiond^. 
When deriving the scaling of the diffusion coefficient with the system size, which involves 
the viscosity of the fluid, Yeh and Hummer investigated the TIP3F® water model com- 
monly used for biomolecular simulations. In the present work, we investigate the diffusion 
coefficient and viscosity of three popular rigid water models of different complexity. The 
SPC/E modeP, involves three interaction sites and is routinely used for the simulation of 
liquid water and electrolyte solutions. The TIP4P/2005 model^, with four interaction sites, 
is usually the best rigid non-polarizable model for the description of the phase diagram 
and the physical properties of bulk solid and liquid water. The Dang-Chang modeP is a 
four-site polarizable model, which allows some transferability from water clusters to liquid 
water and the description of the water-air interface. For these three models, the system 
size-dependence of the diffusion coefficient allows to determine the size-independent value 
and provides an estimate of the viscosity. We also estimate the latter from the Green-Kubo 
relation. 
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II. METHODS 



Model don (A) dou (A) angle (°) eo (kcal/mol) ao (A) qn au (A^) 

SPC/E 1 - 109.471 0.1554 3.1656 0.4238 

TIP4P/2005 0.9572 0.1546 104.52 0.1852 3.1589 0.5564 
Dang-Chang 0.9572 0.215 104.52 0.1825 3.2340 0.5190 1.444 

TABLE I: Geometry and parameters defining the three water models: Bond length don, position of 
the extra M site along the bissector doM, angle between the OH bonds, Lennard- Jones parameters 
on the O atom eo and ao, charge on the H atoms {qo = —2qn for SPC/E, qu = —2qu for the 
other models) and polarizability of the M site. 



The parameters defining the geometry and force field for the three water models are 
summarized in Table |T} The diffusion coefficient are computed using from the mean-squared 
displacement, using the Einstein relation : 

l d(|r(t)-r(0)P) 
^^^^ = /-^^6 d^ 

The "PBC" subscript emphasizes the fact that the use of periodic boundary conditions 

induces a box length dependence of the measured diffusion coefficient ad^ 

2.837A;bT 

DpBc = Do- ^ ; 2 

where rj is the shear viscosity of the solvent. In practice a linear fit of Dpsc" vs 1/L provides 
both Do and an estimate of the viscosity rjpBc- A more usual determination of the viscosity 
consists in using the Green-Kubo relation: 

VGK = / (fTa,3(t)(T«/3(0)) dt (3) 



ksT Jo 

involving the auto-correlation function (ACF) of the off-diagonal components of the stress 
tensor cTo^. 

We performed molecular dynamics simulations in the NVT ensemble, for ranging from 
64 to 4096 for SPC/E, 2048 for TIP4P/2005 and 512 for Dang-Chang (DC). The volume is 
set to ensure a density of 0.998 g.cm~^, and a Nose-Hoover thermostat is used to maintain 
a temperature T = 300 K. Long-range interactions are computed using the Ewald summa- 
tion techniqu^^ and a cut-off is used for short-range interactions. For the DC model, the 



induced dipoles are computed at each time step by minimizing the polarization energy and 
a convergence criterium of 10~^ is used. Simulations of 1 to 10 ns are performed using a 
time step of 1 fs and positions are sampled every 100 fs. The equations of motion for the 
rigid water molecules are integrated using the SHAKE algorithm^^. For simulations with 
the non-polarizable SPC/E and T1P4P/2005 models we use the LAMMP^ code. We also 
compared our results for TIP4P/2005, which contains a massless site, with that obtained us- 
ing Fincham's implicit quaternion algorithmpl instead of SHAKE, and the DLPOLY cod^^. 
Since the same results were found, we only report the ones obtained with SHAKE. Simula- 
tions with the polarizable Dang-Chang model were performed with FIST, the classical part 
of the CP2K simulation package^^. The diffusion coefficients are computed from the slope 
of the mean-square displacement in the 100-300 ps window. Uncertainties on the diffusion 
coefficients for a given system size are estimated using the block averaging methocP^I. The 
uncertainty on the size-independent diffusion coefficient Dq and viscosity estimate rjpBc are 
obtained from the least-square fit of the size- dependent ones to Eq. |2} For the Green-Kubo 
calculations we use the systems with = 256 molecules and collect the components of the 
stress tensor at every time step. The ACF is averaged over the off-diagonal components, and 
the viscosity is computed as the integral over 5 ps (see below). The uncertainty is estimated 
by computing the standard deviation of the values obtained for the different components. 



III. RESULTS 



The results for the diffusion coefficients as a function of system size are reported in 
Figure [1} For the three water models, the scaling predicted by Eq. |2] is observed. The 
corresponding size-independent diffusion coefficients Dq and viscosities rjpBc are summarized 
in Table [IT} All the considered models overestimate the diffusion coefficient compared to the 
experimental result. This is consistent for SPC/E with the result of Kerisit and LiiP^. To 
the best of our knowledge, the size-independent diffusion coefficient of the TIP4P/2005 and 
DC models had not been previously reported. This fact was noticed as a caveat in Ref.'^, 
where it was anticipated that extrapolating to infinite system size would bring the result 
for this model, which is the only rigid one to undersetimate the diffusion coefficient for 
typical box sizes, closer to the experimental one. The extrapolated value Dq turns out to 
be approximately 10% larger. Nevertheless, TIP4P/2005 provides the best agreement with 
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the experimental value. This confirms the conclusion drawn, among rigid non-polarizable 
models, from diffusion coefficients without taking into account the system size-dependencj^. 
In their paper introducing Eq. [2| Yeh and Hummer investigated the TIP3P water modeP and 
obtained Dq = 6.05 10~^ m^s~^, in even worse agreement with the experimental results. The 
polarizable DC model overstimates the diffusion coefficient by less than 20% and performs 
better than SPC/E (30%). 
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FIG. 1: Diffusion coefficient as a function of the inverse box size. The size- independent diffusion 
coefficient Dq is the extrapolation at 1/L — t- 0. The slope provides an estimate of the viscosity 
■qpBC using Eq. [2} 



Model Do (10-9 m^s'^) tjpbc (cP) Vgk (cP) 

SPC/E 2.97±0.05 0.64±0.02 0.68±0.02 

TIP4P/2005 2.49±0.06 0.83±0.07 0.83±0.05 
Dang-Chang 2.72±0.09 0.78±0.06 0.74±0.09 
Exp. 2.3^^ 0.896^ 



TABLE II: Size-independent diffusion coefficient Dq and viscosity {t]pbc from Eq. [2]and r]GK from 
Eq. |3| of the three water models, a: Referencfe^^, b: Referenc^^. 

Before commenting the viscosity estimates rjpBCy we now turn to the Green-Kubo cal- 
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culations. Chen and Smit have recently shown that the viscosity can be computed from 
equihbrium simulations with an accuracy comparable to that obtained with non-equilibrium 
methods, provided that the stress ACF is propertly sampled at short times and that the 
integral Eq. [s] is estimated from times of the order of a few p^. Figure |2] reports the 
stress ACF at very short times and the time-dependent viscosity (running integral of the 
ACF) for the three water models. For SPC/E and TIP4P/2005, the ACF is very similar 
to the one reported by Gonzalez and AbascaP^, with a fast decay and oscillations within a 
few 100 fs which are more pronounced in the TIP4P/2005 case, followed by a slower decay 
which contributes significantly to the integral. The same behaviour is observed for the DC 
water model, for which no such study had been yet reported. As can be seen on Figure |2| 
the choice of an upper limit of 5 ps provides a good estimate of the integral Eq. |3| 

The viscosities tjgk obtained by the Green-Kubo formula and rjpBc obtained from Eq. [2] 
are summarized in Table [Tij For all three models, a good agreement between the two ap- 
proaches is obtained. The values for the SPC/E model are comparable to most of the 
previously reported values, both with equilibrium and non-equilibrium methods: 0.64 cl^, 
0.65 cP'^and 0.67 cP^°. A recent study of water confined in clay nanopores of width larger 
than 4 nm concluded from non-equilibrium molecular dynamics simulations to a viscosity of 
SPC/E water comparable to the bulk value of 0.68-0.70 cP^. Slightly larger values have also 
been reported for bulk SPC/E: 0.72 cF^Sl and 0.73 cF^. The value found for TIP4P/2005 
is slightly smaller than the only reported values in the literature, namely 0.855 and 
0.89 cl^. Nevertheless, we find an excellent agreement between the two methods [rjpBc 
and tjgk)- Since in both references the reported values for SPC/E are also on the larger side 
(0.73 and 0.82 cP, respectively), it is not entirely surprising. To the best of our knowledge, 
no viscosity of the DC model had been previously reported. As for the size-independent dif- 
fusion coefficient, TIP4P/2005 provides the best agreement with the experimental viscosity, 
and the DC model performs better than SPC/E. From the scaling of the diffusion coefficient 
with system size, Yeh and Hummer had found for TIP3P a value of tjpbc = 0.31 ±0.01 cl^, 
in worse agreement with the experimental result, as for the diffusion coefficient. 
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FIG. 2: Stress auto-correlation function and viscosity (see the Green-Kubo equation |3| with (5 
l/ksT) for the three water models. 



IV. CONCLUSION 



We have determined the size-independent diffusion coefficient and shear viscosity of three 
popular rigid water models: SPC/E, TIP4P/2005 and Dang-Chang. A good agreement is 
found for each model between the viscosity determined from the slope of D vs the inverse 
box size and from the Green-Kubo expression. The three considered models overestimate 
the diffusion coefficient and underestimate the viscosity. The TIP4P/2005 model is in better 
agreement with the experimental data for both diffusion and viscosity, followed by the Dang- 
Chang and SPC/E models. 
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This confirms the success of TIP4P/2005 to model not only a large variety of ther- 
modynamic and material properties of bulk water, but also dynamic ones. While it had 
been shown to correctly predict the viscosity of liquid water at ambient conditions, as well 
function of temperature and pressur e * the size-independent diffusion coefficient 
had not been reported previously. Although the prediction of the more complex (and thus 
computationally more costly) Dang-Chang model are in slightly worse agreement with the 
experimental data, one should keep in mind that introducing polarizability should allow for 
a better transferability to other conditions, in particular liquid-solid interfaces and ionic so- 
lutions. The parametrization of a force field compatible with the Dang-Chang water model 
to such situations is currently in progres^. 
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